vonKarmanBeamIRsDoc

Geometrically nonlinear van Karman Beam with internal resonance

Contents

von Karman beam with 1:3 internal resonance

In this example, we consider a clamped-pinned von Karman beam with a support spring at its midspan. The stiffness of the support spring is tund such that 1:3 internal resonance occurs between the first two bending modes. We then extract the forced response curve using SSM reduction.

clear all

Model Setup

nElements = 100;
[M,~,K,fnl,~,Outdof] = build_model(nElements);
outdof = Outdof(2); % the point at mid span
n = length(M);
Kc = K;
Building FE model
Assembling M,C,K matrices
Applying boundary conditions
Solving undamped eigenvalue problem
Getting nonlinearity coefficients
Assembling Tensors
Assembling external force vector

The beam is discretized with 100 elements and then 298 degrees-of-freedom. The dimension of phase space of the full system is 596. With the SSM reduction, the dimension of phase space will be reduced to four. Numerical experiments show that 1:3 near resonance is observed when the stiffness of the support spring is around 37. Next we set the stiffness to be 37 and then perform SSM reduction.

Dynamical System Setup

kLinear = 37;
Kc(outdof,outdof) = K(outdof,outdof)+kLinear;   % linear part
kNonlinear = 0;                                 % nonlinear part
f3_new  = fnl{2};
fnl_new = fnl;
f3_new(outdof,outdof,outdof,outdof) = f3_new(outdof,outdof,outdof,outdof)+kNonlinear; % cubic nonlinerity
fnl_new{2} = f3_new;
C = (2e-4)/9*Kc;

order = 2;
DS = DynamicalSystem(order);
set(DS,'M',M,'C',C,'K',Kc,'fnl',fnl_new);
set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex');

Primiary resonance with IRs

Add forcing - a concentric harmonic force is applied at the midspan

f_0 = zeros(n,1);
f_0(outdof) = 1000;
kappas = [-1; 1];
epsilon   = 0.02;
coeffs = epsilon * [f_0 f_0]/2;
DS.add_forcing(coeffs, kappas);

Create SSM

S = SSM(DS);
set(S.Options, 'reltol', 0.8,'notation','multiindex')

Extract FRC

Here we use the extract_FRC routine to obtain the forced response curve over .

om = eigs(Kc,M,2,'smallestabs');
om = sqrt(om);
freqRange = [0.96 1.05]*om(1);
order = 7;
set(S.FRCOptions, 'nCycle',50000);
set(S.FRCOptions, 'initialSolver','forward');
set(S.contOptions, 'h_max', 0.2, 'PtMX', 1000, 'ItMX',20);
set(S.FRCOptions, 'omegaSampStyle','cocoBD');
set(S.FRCOptions, 'method','continuation ep','outdof',Outdof');
FRC = S.extract_FRC('freq', freqRange, order);
Due to high-dimensionality, we compute only the first 5 eigenvalues with the smallest magnitude. These would also be used to compute the spectral quotients
Assuming a proportional damping hypthesis with symmetric matrices
modal damping ratio for 1 mode is 3.688485e-04
modal damping ratio for 2 mode is 1.106767e-03
modal damping ratio for 3 mode is 2.309722e-03
modal damping ratio for 4 mode is 3.944138e-03
modal damping ratio for 5 mode is 6.019377e-03

 The first 10 nonzero eigenvalues are given as 
   1.0e+02 *

  -0.0001 + 0.3320i
  -0.0001 - 0.3320i
  -0.0011 + 0.9961i
  -0.0011 - 0.9961i
  -0.0048 + 2.0787i
  -0.0048 - 2.0787i
  -0.0140 + 3.5497i
  -0.0140 - 3.5497i
  -0.0326 + 5.4173i
  -0.0326 - 5.4173i

The master subspace has internal resonances: [1  1  3  3]
*****************************************
Calculating FRC using SSM with master subspace: [1  2  3  4]
The master subspace contains the following eigenvalues
 
lambda1 == - 0.0122444 + 33.1964i
 
lambda2 == (-0.0122444) - 33.1964i
 
lambda3 == - 0.110244 + 99.609i
 
lambda4 == (-0.110244) - 99.609i
 
(near) outer resonance detected for the following combinations of master eigenvalues
 They are in resonance with the following eigenvalues of the slave subspace 
 
0*lambda1 + 0*lambda2 + 2*lambda3 + 0*lambda4 == - 0.4801333 + 207.8744i
 
1*lambda1 + 0*lambda2 + 2*lambda3 + 0*lambda4 == - 0.4801333 + 207.8744i
 
0*lambda1 + 0*lambda2 + 3*lambda3 + 1*lambda4 == - 0.4801333 + 207.8744i

..

3*lambda1 + 1*lambda2 + 0*lambda3 + 6*lambda4 == (-3.260961) - 541.7341i
 
sigma_out = 266
(near) inner resonance detected for the following combination of master eigenvalues:
 
0*lambda1 + 2*lambda2 + 1*lambda3 + 0*lambda4 == lambda1
 
1*lambda1 + 0*lambda2 + 1*lambda3 + 1*lambda4 == lambda1
 
2*lambda1 + 1*lambda2 + 0*lambda3 + 0*lambda4 == lambda1

 ..

6*lambda1 + 0*lambda2 + 0*lambda3 + 3*lambda4 == lambda4
 
sigma_in = 266
Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Manifold computation time at order 2 = 00:00:00
Estimated memory usage at order  2 = 9.23E-01 MB
Manifold computation time at order 3 = 00:00:00
Estimated memory usage at order  3 = 2.45E+00 MB
Manifold computation time at order 4 = 00:00:00
Estimated memory usage at order  4 = 3.84E+00 MB
Manifold computation time at order 5 = 00:00:01
Estimated memory usage at order  5 = 9.04E+00 MB
Manifold computation time at order 6 = 00:00:05
Estimated memory usage at order  6 = 1.38E+01 MB
Manifold computation time at order 7 = 00:00:16
Estimated memory usage at order  7 = 2.71E+01 MB

 Run='freqSubint1.ep': Continue equilibria along primary branch.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          4.31e-03  4.95e+01    0.0    0.0    0.0
   1   1  1.00e+00  1.03e-02  8.63e-05  4.95e+01    0.0    0.0    0.0
   2   1  1.00e+00  1.26e-05  1.13e-08  4.95e+01    0.0    0.0    0.0
   3   1  1.00e+00  5.11e-09  3.74e-16  4.95e+01    0.0    0.1    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1         rho2          th1          th2          eps
    0  00:00:00   4.9522e+01      1  EP      3.3196e+01   4.8105e-01   2.1705e-02   6.2609e+00   9.1796e+00   1.0000e+00
   10  00:00:01   4.8851e+01      2          3.2653e+01   3.4855e-01   2.6979e-03   6.2693e+00   9.3302e+00   1.0000e+00
   20  00:00:01   4.8093e+01      3          3.2075e+01   2.4117e-01   5.0683e-04   6.2737e+00   9.3686e+00   1.0000e+00
   24  00:00:01   4.7821e+01      4  EP      3.1869e+01   2.1324e-01   3.0186e-04   6.2748e+00   9.3765e+00   1.0000e+00

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1         rho2          th1          th2          eps
    0  00:00:02   4.9522e+01      5  EP      3.3196e+01   4.8105e-01   2.1705e-02   6.2609e+00   9.1796e+00   1.0000e+00
   10  00:00:02   4.9597e+01      6          3.3428e+01   5.3169e-01   9.4532e-02   6.2079e+00   8.5534e+00   1.0000e+00
   20  00:00:03   4.9267e+01      7  SN      3.3523e+01   5.2095e-01   1.4014e-01   6.1398e+00   7.1645e+00   1.0000e+00
   20  00:00:03   4.9267e+01      8  FP      3.3523e+01   5.2095e-01   1.4012e-01   6.1398e+00   7.1639e+00   1.0000e+00
   20  00:00:03   4.9262e+01      9          3.3523e+01   5.2100e-01   1.3947e-01   6.1410e+00   7.1486e+00   1.0000e+00
   22  00:00:03   4.9213e+01     10  SN      3.3522e+01   5.2356e-01   1.2800e-01   6.1605e+00   6.9623e+00   1.0000e+00
   22  00:00:03   4.9213e+01     11  FP      3.3522e+01   5.2361e-01   1.2786e-01   6.1607e+00   6.9606e+00   1.0000e+00
   30  00:00:04   4.9260e+01     12          3.3650e+01   5.8432e-01   4.8270e-02   6.2464e+00   6.4127e+00   1.0000e+00
   37  00:00:04   4.9500e+01     13  HB      3.3841e+01   6.3662e-01   3.2400e-02   6.2518e+00   6.3167e+00   1.0000e+00
   39  00:00:04   4.9746e+01     14  HB      3.4028e+01   6.8393e-01   2.7200e-02   6.2516e+00   6.2775e+00   1.0000e+00
   40  00:00:04   4.9992e+01     15          3.4211e+01   7.2881e-01   2.4882e-02   6.2503e+00   6.2541e+00   1.0000e+00
   45  00:00:05   5.0867e+01     16  EP      3.4856e+01   8.7616e-01   2.3015e-02   6.2444e+00   6.2074e+00   1.0000e+00
Total time spent on FRC computation upto O(7) = 00:01:15

As an alternative, we can use SSM-ep toolbox to obtain the curve. We also change coordinate representation to Cartesian. The two coordinates yield the same results.

set(S.FRCOptions, 'coordinates', 'cartesian');
set(S.contOptions, 'h_min',1e-2,'h_max',0.1);
S.SSM_isol2ep('isol',[1 2 3 4],order,[1 3],'freq',freqRange,Outdof);
The master subspace contains the following eigenvalues
 
lambda1 == - 0.0122444 + 33.1964i
 
lambda2 == (-0.0122444) - 33.1964i
 
lambda3 == - 0.110244 + 99.609i
 
lambda4 == (-0.110244) - 99.609i
 
(near) outer resonance detected for the following combinations of master eigenvalues
 They are in resonance with the following eigenvalues of the slave subspace 
 
0*lambda1 + 0*lambda2 + 2*lambda3 + 0*lambda4 == - 0.4801333 + 207.8744i
 
1*lambda1 + 0*lambda2 + 2*lambda3 + 0*lambda4 == - 0.4801333 + 207.8744i
 
0*lambda1 + 0*lambda2 + 3*lambda3 + 1*lambda4 == - 0.4801333 + 207.8744i
 
..
 
3*lambda1 + 1*lambda2 + 0*lambda3 + 6*lambda4 == (-3.260961) - 541.7341i
 
sigma_out = 266
(near) inner resonance detected for the following combination of master eigenvalues:
 
0*lambda1 + 2*lambda2 + 1*lambda3 + 0*lambda4 == lambda1
 
1*lambda1 + 0*lambda2 + 1*lambda3 + 1*lambda4 == lambda1
 
2*lambda1 + 1*lambda2 + 0*lambda3 + 0*lambda4 == lambda1
 
0*lambda1 + 2*lambda2 + 2*lambda3 + 1*lambda4 == lambda1

..
 
6*lambda1 + 0*lambda2 + 0*lambda3 + 3*lambda4 == lambda4
 
sigma_in = 266
Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Manifold computation time at order 2 = 00:00:00
Estimated memory usage at order  2 = 9.23E-01 MB
Manifold computation time at order 3 = 00:00:00
Estimated memory usage at order  3 = 2.45E+00 MB
Manifold computation time at order 4 = 00:00:00
Estimated memory usage at order  4 = 3.84E+00 MB
Manifold computation time at order 5 = 00:00:01
Estimated memory usage at order  5 = 9.04E+00 MB
Manifold computation time at order 6 = 00:00:05
Estimated memory usage at order  6 = 1.38E+01 MB
Manifold computation time at order 7 = 00:00:15
Estimated memory usage at order  7 = 2.71E+01 MB

 Run='isol.ep': Continue equilibria along primary branch.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          2.11e-03  4.70e+01    0.0    0.0    0.0
   1   1  1.00e+00  1.85e-03  5.53e-06  4.70e+01    0.0    0.0    0.0
   2   1  1.00e+00  5.43e-06  4.27e-11  4.70e+01    0.0    0.0    0.0
   3   1  1.00e+00  4.54e-11  6.17e-16  4.70e+01    0.0    0.0    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om         Rez1         Rez2         Imz1         Imz2          eps
    0  00:00:00   4.6962e+01      1  EP      3.3196e+01   4.8093e-01  -2.1056e-02  -1.0738e-02   5.2691e-03   1.0000e+00
   10  00:00:00   4.5988e+01      2          3.2509e+01   3.1776e-01  -1.7187e-03  -4.0086e-03   1.4037e-04   1.0000e+00
   20  00:00:00   4.5081e+01      3  EP      3.1869e+01   2.1323e-01  -3.0150e-04  -1.7930e-03   1.4579e-05   1.0000e+00

 STEP      TIME        ||U||  LABEL  TYPE            om         Rez1         Rez2         Imz1         Imz2          eps
    0  00:00:00   4.6962e+01      4  EP      3.3196e+01   4.8093e-01  -2.1056e-02  -1.0738e-02   5.2691e-03   1.0000e+00
   10  00:00:01   4.7313e+01      5          3.3443e+01   5.3055e-01  -5.7972e-02  -4.5720e-02   8.5635e-02   1.0000e+00
   20  00:00:01   4.7394e+01      6          3.3500e+01   5.2364e-01  -2.8209e-03  -7.2212e-02   1.3750e-01   1.0000e+00
   30  00:00:01   4.7423e+01      7          3.3521e+01   5.1583e-01   6.6039e-02  -7.9246e-02   1.2951e-01   1.0000e+00
   35  00:00:02   4.7425e+01      8  SN      3.3523e+01   5.1560e-01   8.9154e-02  -7.4436e-02   1.0812e-01   1.0000e+00
   35  00:00:02   4.7425e+01      9  FP      3.3523e+01   5.1560e-01   8.9155e-02  -7.4436e-02   1.0812e-01   1.0000e+00
   40  00:00:02   4.7425e+01     10  SN      3.3522e+01   5.1963e-01   9.9601e-02  -6.4072e-02   8.0390e-02   1.0000e+00
   40  00:00:02   4.7425e+01     11  FP      3.3522e+01   5.1963e-01   9.9601e-02  -6.4072e-02   8.0388e-02   1.0000e+00
   40  00:00:02   4.7425e+01     12          3.3522e+01   5.2106e-01   9.9953e-02  -6.1415e-02   7.4517e-02   1.0000e+00
   50  00:00:03   4.7458e+01     13          3.3546e+01   5.4515e-01   7.8689e-02  -3.4473e-02   2.6564e-02   1.0000e+00
   60  00:00:03   4.7685e+01     14          3.3706e+01   6.0011e-01   4.1061e-02  -2.0226e-02   3.6559e-03   1.0000e+00
   62  00:00:03   4.7877e+01     15  HB      3.3841e+01   6.3630e-01   3.2382e-02  -1.9990e-02   1.0860e-03   1.0000e+00
   65  00:00:04   4.8142e+01     16  HB      3.4028e+01   6.8359e-01   2.7199e-02  -2.1625e-02  -1.5561e-04   1.0000e+00
   70  00:00:04   4.8659e+01     17          3.4391e+01   7.7114e-01   2.3730e-02  -2.6536e-02  -1.0857e-03   1.0000e+00
   77  00:00:04   4.9320e+01     18  EP      3.4856e+01   8.7550e-01   2.2949e-02  -3.3995e-02  -1.7424e-03   1.0000e+00

FRC in parametrisation space:

FRC in physical space:

Quasi-periodic response

In the above continuation run of fixed points, two Hopf bifurcation points were observed. We switch the continuation of fixed points to the continuation of periodic orbits that are born out of the bifurcation points. Specifically, we use SSM-po toolbox to do that.

HBlab = 15;
set(S.contOptions, 'h_max',0.2,'PtMX',90,'NAdapt',2, 'bi_direct', false,'NSV',2);              % continuation setting
S.SSM_HB2po('po','isol',HBlab,'freq',freqRange,[Outdof; Outdof+n],'saveICs');
 Run='po.po': Continue periodic orbits born from a HB point with label 15 of run isol.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          6.62e-06  4.88e+01    0.0    0.0    0.0
   1   1  1.00e+00  3.71e-06  3.13e-12  4.88e+01    0.0    0.0    0.0
   2   1  1.00e+00  8.72e-11  5.63e-15  4.88e+01    0.0    0.0    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om    po.period          eps
    0  00:00:00   4.8786e+01      1  EP      3.3841e+01   5.8471e+00   1.0000e+00
    2  00:00:00   4.8787e+01      2          3.3842e+01   5.8459e+00   1.0000e+00
..
   90  00:00:22   4.9013e+01     47  EP      3.3907e+01   5.9379e+00   1.0000e+00
Constructing torus in reduced dynamical system
Visualization of torus at (omega,epsilon)=(3.390586e+01,1)

 FRCs from ='po.po': generating torus in physical domain.

Plot of FRC for quasi-periodic orbits

calculate the amplitude of qausi-periodic orbits

runid = 'po';
samps = 1:2:33;
bd   = coco_bd_read([runid,'.po']);
lab  = coco_bd_labs(bd, 'FP');
labs = [1:lab-1 lab+1];
lab  = numel(labs);
x1_tr = zeros(lab,1);
x2_tr = zeros(lab,1);
om_tr = zeros(lab,1);
st    = false(lab,1);
for i=1:lab
    sol = SSM_po_read_solution(runid,labs(i));
    st(i)    = sol.st;
    om_tr(i) = sol.om;
    x1i      = sol.xTr(:,1,:);
    x2i      = sol.xTr(:,2,:);
    x1_tr(i) = norm(x1i(:),'inf');
    x2_tr(i) = norm(x2i(:),'inf');
end

background periodic response

wdir = fullfile(pwd,'data','isol.ep','SSMep.mat');
SSMRes = load(wdir);
FRC = SSMRes.FRC;

plot results of periodic reponse and tori in the same figure

ST = cell(2,1);
ST{1} = {'b--','LineWidth',2}; % unstable
ST{2} = {'b-','LineWidth',2};  % stable
legs = 'SSM-po-unstable';
legu = 'SSM-po-stable';
thm = struct();
thm.SN = {'LineStyle', 'none', 'LineWidth', 2, ...
  'Color', 'cyan', 'Marker', 'o', 'MarkerSize', 8, 'MarkerEdgeColor', ...
  'cyan', 'MarkerFaceColor', 'white'};
thm.HB = {'LineStyle', 'none', 'LineWidth', 2, ...
  'Color', 'black', 'Marker', 's', 'MarkerSize', 8, 'MarkerEdgeColor', ...
  'black', 'MarkerFaceColor', 'white'};

SNidx = FRC.SNidx;
HBidx = FRC.HBidx;
FRC.st = double(FRC.st);
FRC.st(HBidx) = nan;
FRC.st(SNidx) = nan;
figure;
subplot(2,1,1); hold on
plot_stab_lines(FRC.om,FRC.Aout_frc(:,1),FRC.st,ST,legs,legu);
SNfig = plot(FRC.om(SNidx),FRC.Aout_frc(SNidx,1),thm.SN{:});
set(get(get(SNfig,'Annotation'),'LegendInformation'),...
'IconDisplayStyle','off');
HBfig = plot(FRC.om(HBidx),FRC.Aout_frc(HBidx,1),thm.HB{:});
set(get(get(HBfig,'Annotation'),'LegendInformation'),...
'IconDisplayStyle','off');
ylabel('$||w_{0.25l}||_{\infty}$','Interpreter','latex');
set(gca,'FontSize',14);
grid on, axis tight;
legend boxoff;
subplot(2,1,2); hold on
plot_stab_lines(FRC.om,FRC.Aout_frc(:,2),FRC.st,ST,legs,legu);
SNfig = plot(FRC.om(SNidx),FRC.Aout_frc(SNidx,2),thm.SN{:});
set(get(get(SNfig,'Annotation'),'LegendInformation'),...
'IconDisplayStyle','off');
HBfig = plot(FRC.om(HBidx),FRC.Aout_frc(HBidx,2),thm.HB{:});
set(get(get(HBfig,'Annotation'),'LegendInformation'),...
'IconDisplayStyle','off');
xlabel('$\Omega$','Interpreter','latex');
ylabel('$||w_{0.5l}||_{\infty}$','Interpreter','latex');
grid on;

figure(gcf);
subplot(2,1,1);
plot(om_tr(st), x1_tr(st), 'o','MarkerSize', 8, 'MarkerEdgeColor', ...
  'blue', 'MarkerFaceColor', 'blue','DisplayName','SSM-tor-stable');
plot(om_tr(~st), x1_tr(~st), 'o','MarkerSize', 8, 'MarkerEdgeColor', ...
  'red', 'MarkerFaceColor', 'red','DisplayName','SSM-tor-unstable');
lowb = 0.9*min(FRC.Aout_frc(HBidx,1));
upb  = 1.05*max(x1_tr);
axis([min(FRC.om(HBidx))-0.05,max(FRC.om(HBidx))+0.05,lowb,upb]);
subplot(2,1,2);
plot(om_tr(st), x2_tr(st), 'o','MarkerSize', 8, 'MarkerEdgeColor', ...
  'blue', 'MarkerFaceColor', 'blue','DisplayName','SSM-tor-stable');
plot(om_tr(~st), x2_tr(~st), 'o','MarkerSize', 8, 'MarkerEdgeColor', ...
  'red', 'MarkerFaceColor', 'red','DisplayName','SSM-tor-unstable');
lowb = 0.9*min(FRC.Aout_frc(HBidx,2));
upb  = 1.05*max(x2_tr);
axis([min(FRC.om(HBidx))-0.05,max(FRC.om(HBidx))+0.05,lowb,upb]);
set(gca,'FontSize',14);
legend boxoff;